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Abstract The Helioseismic and Magnetic Imager onboard the Solar Dynamics Observa- 
tory (SDO/HMI) provides continuous full-disk observations of solar oscillations. We de- 
velop a data-analysis pipeline based on the time-distance helioseismology method to mea- 
sure acoustic travel times using HMI Doppler-shift observations, and infer solar interior 
properties by inverting these measurements. The pipeline is used for routine production of 
near-real-time full-disk maps of subsurface wave-speed perturbations and horizontal flow 
velocities for depths ranging from 0 to 20 Mm, every eight hours. In addition, Carring- 
ton synoptic maps for the subsurface properties are made from these full-disk maps. The 
pipeline can also be used for selected target areas and time periods. We explain details of 
the pipeline organization and procedures, including processing of the HMI Doppler obser- 
vations, measurements of the travel times, inversions, and constructions of the full-disk and 
synoptic maps. Some initial results from the pipeline, including full-disk flow maps, sunspot 
subsurface flow fields, and the interior rotation and meridional flow speeds, are presented. 

Keywords Sun: helioseismology • Sun: oscillations • Sun: SDO 


1. Introduction 

The Helioseismic and Magnetic Imager onboard the Solar Dynamics Observatory (SDO/ 
HMI: Schou et al, 201 1) observes the solar full-disk intensity, Doppler velocity, and vector 
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magnetic field of the photosphere with high spatial resolution and high temporal cadence. 
Similar to the Michelson Doppler Imager (MDI: Scherrer et al., 1995), an instrument on- 
board the Solar and Heliospheric Observatory (SOHO), the HMI Dopplergrams are pri- 
marily used for helioseismic analysis to investigate the interior structure and dynamics of 
the Sun. Helioseismology data-analysis pipelines are planned for near-real-time analyses of 
the observations in order to provide the analysis results to the helioseismology and solar 
physics communities. The time-distance analysis pipeline is one of the pipelines for local 
helioseismology studies, and other pipelines include ring-diagram analysis and far-side ac- 
tive region imaging. The time-distance pipeline is designed for the routine production of 
nearly full-disk subsurface wave-speed perturbations and horizontal flow fields every eight 
hours, as well as synoptic flow maps for every Carrington rotation. It can also be used to 
analyze specific target areas and time periods. 

Time-distance helioseismology was first introduced by Duvall et al. (1993, 1996), and it 
has developed rapidly since then. Different inversion techniques were introduced and tested. 
The LSQR algorithm, introduced by Kosovichev (1996) and used later by Zhao, Kosovichev, 
and Duvall (2001), solves the inversion problem in the least-squares sense in the spatial 
domain by an iterative approach. The Multi-Channel Deconvolution (MCD) method, intro- 
duced by Jacobsen et al. (1999) and widely used in later studies ( e.g . Couvidat et al., 2004), 
solves the least-squares problems in the Fourier domain. Later, Couvidat et al. (2005) ap- 
plied a horizontal-regularization procedure for this inversion technique. More recently, an 
optimally localized averaging (OLA) inversion scheme was introduced to study the solar 
subsurface flow fields (Jackiewicz, Gizon, and Birch, 2008). 

Different types of sensitivity kernels, which describe the relationship between the travel 
times and interior properties, were also introduced and used in the time-distance inversion 
problems. Kosovichev (1996) first used ray-path approximation kernels, Jensen, Jacobsen, 
and Christensen-Dalsgaard (2000) introduced Fresnel-zone kernels; and Birch and Koso- 
vichev (2000), Birch, Kosovichev, and Duvall (2004), and Birch and Gizon (2007) investi- 
gated Born- approximation kernels for both sound-speed structures and flow fields. Couvidat, 
Birch, and Kosovichev (2006) compared subsurface sound-speed perturbation structures in- 
ferred from these different types of kernels, and found that the inversion results obtained 
with the different kernels were basically consistent. 

Important results on the solar interior properties have been obtained from the time- 
distance studies as well as from other local helioseismology techniques (e.g., Komm et 
al, 2004; Lindsey and Braun, 2000). The introduction that follows is limited to only 
time-distance results due to the scope of this paper. On global scales, poleward merid- 
ional flows were found below the photosphere (Giles et al., 1997), and solar-cycle depen- 
dent meridional flow variations were also investigated and discussed (Chou and Dai, 2001; 
Beck, Gizon, and Duvall, 2002; Zhao and Kosovichev, 2004). On local scales, subsurface 
sound-speed perturbations and flow fields were derived for supergranulation (Kosovichev 
and Duvall, 1997; Duvall et al., 1997; Duvall and Gizon, 2000; Sekii et al., 2007; Jack- 
iewicz, Gizon, and Birch, 2008) and for sunspots (Kosovichev, Duvall, and Scherrer, 2000; 
Gizon, Duvall, and Larsen, 2000; Zhao, Kosovichev, and Duvall, 2001; Couvidat, Birch, 
and Kosovichev, 2006; Zhao, Kosovichev, and Sekii, 2010). Additionally, time-distance 
helioseismology was used to detect the emergence of active regions before their appear- 
ances in the photosphere (Kosovichev, Duvall, and Scherrer, 2000; Jensen et al., 2001; 
Zharkov and Thompson, 2008), to image large active regions on the far side of the Sun 
(Zhao, 2007; Ilonidis, Zhao, and Hartlep, 2009), and to measure sound-speed perturba- 
tions in the tachocline (Zhao et al., 2009). These results are important for space-weather 
forecasting and understanding the mechanisms for the generation of solar magnetism. The 
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time-distance helioseismology pipeline analyses, based on the high spatial-resolution and 
high temporal-cadence observations from HMI, will greatly advance our knowledge of the 
interior processes and their connections with solar activity above the photosphere. 

However, one should keep in mind that the physics of solar oscillations in the turbulent 
magnetized plasma is very complicated, and that the helioseismic techniques are still in the 
process of being developed. Because of limited knowledge of the wave physics and com- 
plexity of the MHD turbulence, there may be systematic uncertainties in the local helioseis- 
mology inferences, particularly in strong magnetic-field regions of sunspots. For example, 
Lindsey and Braun (2005 a, 2005b) argued that the outgoing and ingoing travel-time asym- 
metries observed in sunspot areas might be caused by a “shower-glass effect”. Schunker et 
al. (2005) found that the inclined magnetic field in sunspot penumbrae might cause varia- 
tions of measured acoustic travel times. Zhao and Kosovichev (2006) found that this inclined 
magnetic-field dependent effect does not exist in the measurements obtained from the MDI 
intensity observations. Later, Rajaguru et al. (2006) found that the non-uniform acoustic 
power distribution in the photosphere also contributed to measured travel-time shifts in ac- 
tive regions if a phase-speed filtering procedure was applied. This effect was then studied by 
Parchevsky, Zhao, and Kosovichev (2008) and Hanasoge et al. (2008) numerically, and also 
by Nigam and Kosovichev (2010) analytically. More recently, Gizon et al. (2009) derived 
a sunspot’s subsurface flow fields after applying ridge filtering, and their inferred results 
did not agree with the previously inverted results with the use of phase-speed filtering. Us- 
ing high spatial-resolution observations from Hinode, Zhao, Kosovichev, and Sekii (2010) 
showed that the principal results on sunspot structure did not depend on the use of phase- 
speed filtering. However, significant systematic uncertainties in sunspot seismology remain 
and need to be understood, and these are being actively studied by the use of numerical sim- 
ulations. For a recent review of the sunspot seismology and uncertainties see Kosovichev 
( 2010 ). 

Despite the ongoing discussions of accuracy of time-distance measurements and inter- 
pretation of inversion results, it is useful to provide the measured travel times and the inver- 
sion results for investigations of structures and flows below the visible surface of the Sun. As 
the flow chart of the pipeline in Figure 1 shows, we apply phase-speed filtering to the HMI 
data, compute cross-covariances of the oscillations, and use two different travel-time fitting 
procedures to derive the acoustic travel times. We then perform inversions using two differ- 
ent sets of travel-time sensitivity kernels, based on the ray-path and Born approximations, to 
infer subsurface wave-speed perturbations and flow fields, using the MCD inversion method. 
We provide online access to the measured travel times, full-disk subsurface wave-speed per- 
turbations and flow maps calculated every eight hours. In addition, we provide synoptic 
flow maps for each Carrington rotation. In this article, we describe details of the acoustic 
travel-time measurement procedure in Section 2, and the inversion procedure in Section 3. 
We describe the pipeline data products and present initial HMI results in Section 4, and 
summarize our work in Section 5. 


2. Acoustic Travel-Time Measurement 

2.1. Tracking and Remapping 

SDO/HMI continuously observes the full-disk Sun, providing Doppler velocity, continuum 
intensity, line-depth, line-width, and magnetic-field maps with a 45-second cadence, and 
also vector magnetic-field measurements with a cadence of 12 minutes. Each full-disk image 
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Figure 1 Flow chart for the HMI time-distance helioseismology data-analysis pipeline. 


has 4096 x 4096 pixels with a spatial resolution of 0.504 arcsec pixel - 1 ( i.e ., approximately, 
0.03 heliographic degree pixel - 1 at the solar disk center). The Doppler observations are 
primarily used for helioseismic studies. The observing sequences, algorithms for deriving 
Doppler velocity and magnetic field, and other instrument calibration issues, are discussed 
by Schou et al. (2011). 

As illustrated in Figure 1, the primary input for the pipeline is Dopplergrams, although 
in principle, the HMI intensitygrams and line-depth data can also be analyzed in the same 
manner. Users of the pipeline can select specific areas for analysis, preferably within 60° 
from the solar disk center. In practice, the users provide the Carrington longitude and latitude 
of the center of the selected area, and the mid-time of the selected time period, then the 
pipeline code selects an area of roughly 30° x 30° centered at the given coordinate, and 
for a time interval of eight hours with the given time as the middle point. The data for this 
selected area and time period are then tracked to remove solar rotation, and remapped into 
the heliographic coordinates using Postel’s projection (also known as azimuthal equidistant 
projection) relative to the given area center. Normally, the tracked area consists of 5 12 x 5 12 
pixels with a spatial sampling of 0.06° pixel -1 ; and the temporal sampling is the same 
as the observational cadence. Cubic interpolation is used for the pixels not located on the 
observational grid. 

Figure 2 shows typical HMI k-to and time-distance diagrams obtained from one tracked 
and remapped area. The selected area covers 30° in latitude and has an apparent differential 
rotation over this span. However, for fast computations, only one uniform tracking rate cor- 
responding to the Snodgrass rotation rate at the center of this area is used. The Snodgrass 
rotation rate is 2.851 — 0.343 sin 2 cj) — 0.474 sin 4 4> prads -1 , where </> is latitude (Snodgrass 
and Ulrich, 1990). The uniform tracking rate of the selected area results in a differential 
rotation velocity in the inverted horizontal velocity fields. This differential rotation velocity 
is removed from the full-disk flow maps after averaging over the whole Carrington rotation, 
and only the residual flow fields are given as the final results (see Section 4.1). But for the 
user-selected areas, the differential rotation is kept in the inversion results. 
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Figure 2 Typical power 
spectrum ( k-a > ) diagram (upper) 
and time-distance 
(cross-covariance) diagram 
(lower) made from eight hours of 
HMI Doppler observations. 
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2.2. Computing Cross-Correlations and Fitting for Travel Times 

Each tracked and remapped Dopplergram datacube is filtered in the 3D Fourier domain. 
Solar convection and /-mode oscillation signals are removed first, and then phase-speed 
filtering is applied following the procedures prescribed by Couvidat et al. (2005). For the 
travel-time measurements, for each central point we select 1 1 annuli with various radii and 
widths chosen from our past experience with MDI analyses. All of the phase-speed filter- 
ing parameters, including the central phase speed, filter width, and the corresponding inner 
and outer annulus radii are presented in Table 1. The phase-speed filter is a Gaussian func- 
tion of the wave’s horizontal phase speed. It selects wave packets traveling between the 
central points and the annuli for the selected distances. This filter helps improve the signal- 
to-noise ratio, and also it removes leakage from low-degree oscillations. After the filtering, 
the data are transformed back to the space-time domain for cross-covariance computations. 
Two different fitting methods are used to derive the acoustic travel times from the cross- 
covariances: a Gabor-wavelet fitting (Kosovichev and Duvall, 1997), and a cross-correlation 
method based on seismology algorithms (Zhao and Jordan, 1998) adopted by Gizon and 
Birch (2002; GB algorithm hereafter). A detailed description of the filtering procedure, the 
cross-covariance computations, the two types of travel-time fittings, comparison of the travel 
times derived from the two fitting methods, and the measurement error estimates is given 
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Table 1 Phase-speed filtering 

parameters used for the selected Annulus No. Annulus range Phase speed FWHM 

travel distances (annulus ranges). (heliographic degree) (pHz/£) (pHz/£) 


1 

0.54-0.78 

3.40 

1.0 

2 

0.78-1.02 

4.00 

1.0 

3 

1.08-1.32 

4.90 

1.25 

4 

I.44-I.8O 

6.592 

2.149 

5 

1.92-2.40 

8.342 

1.351 

6 

2.40-2.88 

9.288 

1.183 

7 

3.12-3.84 

10.822 

1.895 

8 

4.08-4.80 

12.792 

2.046 

9 

5.04-6.00 

14.852 

2.075 

10 

6.24-7.68 

17.002 

2.223 

11 

7.68-9.12 

19.133 

2.039 


by Couvidat et al. (2010). In particular, it has been shown there that the two definitions 
of the acoustic travel times and the two fitting methods give generally consistent results in 
the quiet-Sun regions, but that they may give different results in active regions. The differ- 
ences are currently not well understood, but in the pipeline we implement both travel-time 
definitions. 

In the Postel-projected maps, the exact distance between two arbitrary points cannot 
be calculated using the Cartesian coordinates. Thus, when an annulus is selected around a 
given location, some additional computations are needed to determine the exact great-circle 
distance between the points. The formula to determine the great-circle distance is 

A = arccos(sin0i sin02 + cos^i cos02cos(0i —((>2)), (1) 

where (6fi, 0i) and (62, 4> 2 ) are the heliospheric longitude and latitude coordinates for the 
two separate locations. 

To facilitate the inversions for subsurface flow fields, each annulus is divided into four 
quadrants representing the North, South, East, and West directions (Kosovichev and Duvall, 
1997). So, for each annulus and each travel-time fitting method, we obtain six measurements 
of acoustic travel times, corresponding to the outgoing and ingoing, East-, West-, South-, 
and North-going directions. We then combine these travel times to obtain one map for the 
mean travel time and three maps for the travel-time differences. These travel times are: r mea n 
(average of outgoing and ingoing travel times), r 0 , (difference of outgoing and ingoing travel 
times), r we (difference of West- and East-going travel times), and r ns (difference of North- 
and South-going travel times). These four travel-time maps for each annulus are archived 
and available through the HMI Data Record Management System (DRMS). 


3. Subsurface Wave-Speed Perturbation and Flow-Field Inversions 

The acoustic travel times are derived by two different fitting methods: the Gabor- wavelet 
function and the GB algorithm. Then, as illustrated in Figure 1 , to infer the subsurface wave- 
speed perturbations and flow velocities, the Gabor-wavelet travel times are inverted using 
the ray-path approximation sensitivity kernels, and the GB travel times are inverted using 
the Born- approximation sensitivity kernels. The Born- approximation kernels are calculated 
based with the filter and window parameters of the GB fits. 
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3.1. Inversions 

Both the ray-path and Born- approximation kernels have been used in previous time-distance 
studies (see, e.g., Zhao, Kosovichev, and Duvall, 2001; Couvidat, Birch, and Kosovichev, 
2006). Details of the kernel calculations and their comparisons will be given in a separate 
article. 

We employ the MCD inversion method (Jacobsen et ah, 1999) with a horizontal regu- 
larization (Couvidat et ai, 2005). For the wave-speed perturbation inversions, the linearized 
equation relating the measured mean travel times and the subsurface wave-speed perturba- 
tions are 

(2) 

ijk 

where 8s ljk is the relative wave-speed perturbation 8c i j k /c i j k approximated by piece-wise 
constant functions on the inversion grid, and A^ k v is a matrix of the discretized sensitivity 
kernel. Here, k and /x label the coordinates of the central points of the annuli in the observed 
areas, v labels the different annuli, and i, j, and k are the indices for the discretized three- 
dimensional space for inversions. Usually, the horizontal coordinates of the inversion grid 
(i and j indices) are selected at the same locations as the central travel-time measurement 
points (A. and /x indices). In the first-order approximation, the sensitivity kernels are calcu- 
lated for a spherically symmetric solar model and do not depend on the position on the solar 
surface. Therefore, in this case Equation (2) is actually equivalent to a convolution in the 
horizontal domain, which can be simplified as a direct multiplication in the Fourier domain: 

<5f l '(/cy,/qJ = 'Y^ J A v k {K k ,ic ll )8s k {Kx,K IJ ,), (3) 

k 

where St , A, and 8s are the 2D Fourier transforms of St, A, and 8s, respectively; k is the 
same as in Equation (2); k,_ and K k are the wavenumbers in the Fourier domain correspond- 
ing to k and fx of the spatial domain. For each (k-,. k^), the equation in the Fourier domain 
is a matrix multiplication: 

d = Gm , (4) 

where 

d = {8 t v {kx,k ij )}, G={A v k {Kx,K tl )), m = {8s k (K),,K^)}- 

Thus, we have a large number of linear equations describing the depth dependence of the 
Fourier components, and each linear equation can be solved in the least-squares sense. After 
all these equations are solved, and m is obtained for each (k> , ) , the values of 8s ijk are 
calculated by the inverse 2D Fourier transform. 

Equation (4) is ill-posed, and regularization is required to obtain a smooth solution. The 
regularized least-squares algorithm is formulated as 

min{ | (z/ — Gm) + k 2 (K)\\Lm\\\}, (5) 

where || . . . || 2 denotes the F2-norm, L is a regularization operator, and a (at) is a regulariza- 
tion parameter. We choose L to be a diagonal matrix whose elements are the inverse of the 
square root of the spatial sampling Az at each depth. Such weighting is necessary because 
the grid in the vertical direction is chosen to be approximately uniform in the acoustic depth, 
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meaning that the spatial sampling of deep layers is larger than the sampling of shallower lay- 
ers. The regularization parameter is taken in the form of X 2 (k) = a 2 + X\(k), where X v and 
A/ ? are vertical and horizontal regularization parameters. The purpose of introducing A/, is 
to damp the high-wavenumber components that may lead to noise amplification. Following 
the discussion of Couvidat et al. (2005), we choose A./, = X 2 K 2 with A 2 as a constant. 

Because the regularization is applied in the Fourier domain, it is quite difficult to use dif- 
ferent regularization parameters for different horizontal locations of the same region. Some- 
times different regularization parameters are needed, because in active regions the noise 
level may be different from the quiet-Sun regions, as we discuss in Section 3.3. Thus, it is 
necessary to implement into the analysis pipeline another inversion technique, the LSQR 
algorithm, which solves Equation (2) in the space domain by an iterative approach. This 
implementation is currently under development. 

3.2. Inversion Depth and Validation of Inversions 

For both the wave-speed and flow-field inversions, and for both the ray-path and Born- 
approximation inversions, we select a total of 1 1 inversion depths as follows: 0-1, 1—3, 3 — 
5, 5-7,7-10, 10-13, 13-17, 17-21, 21-26, 26-30, and 30-35 Mm. There are a total 
of 1 1 depth intervals. The inversion results provide the wave-speed perturbations and flow 
velocities averaged in these layers. Due to the lack of acoustic-wave coverage in the deep 
interior, the reliability of inversion results decreases with the depth. Thus, only inversion 
results for the depths shallower than 20 Mm are included in the pipeline output. This may 
change in the future when more confidence is gained in the deeper interior inversion results. 

In recent years, several studies have been carried out to validate the time-distance mea- 
surements and inversions. To validate the derived subsurface flow fields, Georgobiani et 
al. (2007) and Zhao et al. (2007) have analyzed realistic solar-convection simulations and 
found satisfactory inversion results for the shallow layers covered by the simulations. Vali- 
dations of the wave-speed perturbation inversions based on numerical simulations with pre- 
set structures have also been performed. Meanwhile, numerical simulations for magnetic 
structures with flows are also under development (Rempel, Schiissler, and Knolker, 2009; 
Stein et al., 2011; Kitiashvili et al., 2011). Validations of the time-distance helioseismology 
techniques will be carried out as well using these simulations. 

Cross-comparisons between different helioseismology techniques, e.g. comparing the 
mean rotation speed from our pipeline analysis with global helioseismology results, and 
comparing the subsurface flow fields with results from ring-diagram analyses, will also be 
important for the validation. 

3.3. Error Estimate 

There are two types of errors in the pipeline results: systematic errors due to our limited 
knowledge of the wave physics in the magnetized turbulent medium and simplified mathe- 
matical formulations, and statistical errors, which are mostly due to the stochastic nature of 
the solar oscillations (“realization noise”). Flere, we only discuss the statistical errors. 

To estimate the errors in the inversion results, we need first to estimate the uncertain- 
ties in the measured acoustic travel times. Flere, we estimate the measurement uncertainties 
following the prescriptions of Jensen, Duvall, and Jacobsen (2003) and Couvidat, Birch, 
and Kosovichev (2006). We select 25 quiet-Sun regions, and use the rms variation of mean 
travel times for different measurement distances as an error estimate for the travel times. The 
estimated uncertainties for the Gabor-wavelet fitting are given in Table 2, and the uncertain- 
ties obtained for the GB algorithm are similar, but slightly larger for short distances and 
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Table 2 Measured mean travel times and uncertainties for both quiet-Sun regions and active regions. 


Annulus No. 

Mean travel time 
(min) 

Uncertainty for quiet regions 
(min) 

Uncertainty for active regions 
(min) 

1 

11.87 

0.062 

0.17 

2 

18.82 

0.061 

0.25 

3 

21.76 

0.11 

0.26 

4 

25.85 

0.11 

0.19 

5 

28.69 

0.14 

0.15 

6 

31.18 

0.14 

0.14 

7 

35.07 

0.10 

0.10 

8 

38.86 

0.12 

0.11 

9 

42.46 

0.11 

0.093 

10 

46.63 

0.14 

0.11 

11 

50.26 

0.15 

0.11 


slightly smaller for long distances. Active regions have different measurement uncertainties. 
To estimate these, we selected a relatively stable sunspot, NOAA AR 1 1092, from 2 through 
5 August 2010, and we assumed that the sunspot did not change during this period. Then 
we use the rms of the travel times measured inside the sunspot as an error estimate. Due to 
the evolution of the sunspot, this approach overestimates the measurement uncertainties, but 
can still give us an approximate estimate of measurement errors. These error estimates are 
presented in Table 2 as well. 

Inversions are then performed for the same quiet-Sun regions and the active region to 
estimate the statistical errors in inversion results. Following the same approach, the rms of 
inverted wave-speed perturbations is assumed as the statistical error. The error estimates for 
both the quiet-Sun and active regions are shown in Table 3. Because supergranular flows are 
dominant in the flow fields, it is difficult to estimate errors of the inverted velocity for the 
quiet Sun by this approach. Instead, we estimate the velocity errors based on the rotational 
velocity profile, which has little change in a time scale of a few days. While the errors for 
the wave-speed perturbation inside active regions are roughly twice those for the quiet-Sun 
regions, the velocity errors inside active regions sometimes are seven times larger than the 
errors in the quiet Sun. 


4. Data Products and Initial Results from HMI 

The time-distance data-analysis pipeline is used for the routine production of nearly real- 
time full-disk (actually, nearly full-disk covering a 120° x 120° area on the solar disk) wave- 
speed perturbation and flow field maps every eight hours. These maps are then used to 
construct the corresponding synoptic maps for each Carrington rotation. The pipeline can 
also be used for specific target areas, such as active regions. In this section, we introduce the 
data products from this pipeline and some initial results from it. 

4.1. Routine Production: Full-Disk and Synoptic Maps 

For each day of HMI observations, we select three eight-hour periods: 00:00-07:59 UT, 
08:00-15:59 UT, and 16:00-23:59 UT. For each analysis period, we select 25 regions, 
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Table 3 Error estimates for the relative wave-speed perturbation and horizontal-velocity inferences in the 
quiet-Sun (QS) and active regions (AR). 


Depth 

(Mm) 

Wave speed for QS 

Velocity for QS 
(ms -1 ) 

Wave speed for AR 

Velocity for AR 
(ms -1 ) 

0-1 

2.8 x 10 -3 

7.8 

6.3 xlO -3 

58.3 

1-3 

4.1 x 10 -3 

7.5 

10.9 x 10 -3 

56.4 

3-5 

6.4 x 10 -3 

8.9 

8.7 x 10 -3 

51.1 

5-7 

4.6 x 10 -3 

9.4 

9.7 x 10 -3 

45.1 

7-10 

4.7 x 10 -3 

13.1 

6.7 x 10 -3 

34.5 

10-13 

3.7 x 10 -3 

12.9 

3.1 xlO -3 

28.1 


Figure 3 Schematic plot 
showing how areas are selected 
for a routine calculations of the 
full-disk wave-speed and flow 
maps. Not all of the 25 selected 
areas are shown. 


90°N 



90°S 


with the central locations at 0°, ±24°, and ±48° in both longitude and latitude, where the 
longitude is relative to the central meridian at the mid-time of the selected period. Figure 3 
shows locations of these areas on the solar disk. The total number of areas is 25: five rows 
and five columns. Due to the Postel’s projection, the boundaries of these areas are often not 
parallel to the latitude or longitude lines. It is also evident that many areas overlap, some 
areas overlap twice, and some overlap four times. The travel times and inversion results are 
averaged in these overlapped areas. However, in the areas close to the solar limb, the fore- 
shortening effect may become non-negligible, but the role of this is not yet systematically 
studied. Users of these maps are urged to be cautious when using the pipeline results in the 
areas close to the limb. 

For each full-disk map and each synoptic map, the East-West velocity | v x | , the North- 
South velocity [u y ], and wave-speed perturbation [<5c/c] in each depth layer are derived 
with a horizontal spatial sampling of 0.12° pixel -1 . For each of the 25 areas, the inversion 
results are first obtained in the Postel-projection coordinates, and then converted into the 
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Figure 4 A map of horizontal flow divergence for a depth range of 1 - 3 Mm and a time period of 
00:00-07:59 UT 19 May 2010. The display scale is from —6.2 x 10“ 4 to 6.2 x 1CP 4 s^ 1 . White areas 
with positive divergence correspond to supergranulation. Note that supergranules appear larger at high lati- 
tudes because of the rectangular longitude-latitude map projection. 


longitude-latitude coordinates. This coordinate conversion is basically the inverse procedure 
of transforming the observed data into the Postel-projection coordinates for the travel-time 
measurements. Cubic spline interpolation is employed. The results in high-latitude regions 
are oversampled. After the coordinate transformation, the overlapped areas are averaged, 
and the final statistical errors are estimated for the whole procedure starting from the travel- 
time measurements. This includes all potential errors from the interpolation and remapping. 
The final full-disk results are saved on a uniform longitude-latitude coordinate grid, so each 
horizontal image of the subsurface layers has a total of 1000 x 1000 pixels covering 120° in 
both longitudinal and latitudinal directions. This coordinate choice is convenient for merging 
and analyzing results, but unavoidably distorts maps in high latitude areas. 

Figure 4 shows an example of a full-disk map of the subsurface horizontal flow diver- 
gence calculated from v x and v y at the depth range of 1 - 3 Mm. The positive-divergence 
areas correspond to supergranulation. Such maps at various depths with continuous temporal 
coverage can be valuable for studying the structure and evolution of the supergranulation. 
Figure 5 displays the subsurface horizontal flow fields with full spatial resolution for a region 
located at the center of the map in Figure 4. Supergranular flows can be easily identified. 
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Figure 5 A sample of subsurface horizontal flow fields with full spatial resolution at the depth of 0 - 1 Mm. 
This area is sampled at the center of the map shown in Figure 4. The background image shows the line-of-sight 
magnetic field measured by HMI, with red as positive and blue as negative polarities. The range of the 
magnetic field is from —80 to 80 Gauss. 


From the full-disk wave-speed perturbation and flow maps obtained every eight hours, 
we select an area 13.2° wide in longitude, i.e. —6.6° to 6.6° from the central meridian, to 
construct the synoptic wave-speed perturbation and flow maps. Since the Carrington rotation 
rate corresponds to a shift of approximately 4.4° every eight hours, each location in the con- 
structed synoptic maps is averaged roughly three times (i.e. one whole day). The resultant 
synoptic map for each depth consists of 3000 x 1000 pixels. Since such a map is difficult 
to display, we show in Figure 6 a binned-down synoptic large-scale flow map obtained for 
the depth of 1 - 3 Mm for Carrington Rotation 2097 during the period from 20 May to 16 
June, 2010. The vectors in the figure are obtained by averaging the flow fields in areas of 
15° x 15° with a 3° sampling rate. The maps of this type are similar to the subsurface flow 
maps obtained from the ring-diagram analysis (Flaber et al, 2002). From Figure 6, it can be 
found that the large-scale flows often converge around magnetic regions. 

Normally, the full-disk flow maps and the synoptic flow maps are provided as residual 
flow maps after subtracting the flows averaged over the entire Carrington rotation. The sub- 
tracted average maps contain the differential rotation, meridional flows, and possibly some 
systematic effects. Figure 7 shows examples of the subsurface differential rotation speed and 
meridional flow speed obtained by averaging the calculations for Carrington Rotation 2097. 
These results are also provided online together with the synoptic flow maps. 
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Figure 6 Synoptic map of large-scale horizontal flows at the depth of 1 - 3 Mm for CR 2097. The back- 
ground image is the corresponding line-of-sight magnetic field, with red as positive and blue as negative 
polarities, in the range of —50 to 50 Gauss. 


Figure 7 Averaged rotation 
(upper) and meridional flow 
speeds (lower) at selected depths 
for Carrington Rotation 2097. 
The rotation speed is relative to 
the constant Carrington rotation 
rate. 




4.2. Target Areas 

As already mentioned in Section 2, the pipeline can also be run for specific target areas and 
specific time intervals. The pipeline users are required to provide the Carrington coordinates 
of the center of the target area, and the mid-time of the time interval. 
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Figure 8 Subsurface flow fields of AR 11072 at the depth of 1 -3 Mm, obtained from 16:00-23:59 UT 23 
May 2010. The background image is line-of-sight magnetic field, with red positive and blue negative. The 
image is displayed with a scale from —1000 to 1000 Gs. The arrows are displayed after a 2 x 2 rebinning, 
and the longest arrow represents a speed of 300 ms -1 . 


Figure 8 shows an example of a small part (approximately 1 /9) of a target area, with an 
active region, AR 11072, included in it. The subsurface flow field, at the depth of 1 - 3 Mm, 
is displayed after a 2 x 2 rebinning. The displayed time period, 16:00-23:59 UT on 23 May 
2010, is approximately 2.5 days after the start of emergence of this active region that was 
still growing. Our results clearly show subsurface outflows around the leading sunspot, and 
some converging flows inside it. Comparing with the previous results of Kosovichev (2009) 
and Zhao, Kosovichev, and Sekii (2010), one may conclude that the subsurface flow fields of 
active regions evolve with the evolution of active regions. There may be prominent outflows 
around sunspots during their early growing phase and also their decaying phase, but there 
may also be strong converging flows during the stable period. Our continuous monitoring 
of the solar full-disk subsurface flows may give us an opportunity to statistically study the 
changes of the sunspot’s subsurface flows with the sunspot evolution. 


5. Summary 

We have developed a time-distance helioseismology data-analysis pipeline for SDO/HMI 
Doppler observations. This pipeline performs acoustic travel-time measurements based on 
two different methods, and conducts inversions based on two different sensitivity kernels 
calculated in the ray-path and Bom approximations. The pipeline gives nearly real-time 
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routine products of full-disk wave-speed perturbations and flow-field maps in the range of 
depth 0-20 Mm every eight hours, and provides the corresponding synoptic wave-speed 
perturbation and flow-field maps for each Carrington rotation. The averaged rotation speed 
and meridional flow speed are also provided separately for each rotation. In addition to these 
routine production, the pipeline can also be used for analysis of specific target areas for 
specific time intervals. This data-analysis pipeline will provide important information about 
the subsurface structure and dynamics on both local and global scales, and its continuous 
coverage through years to come will be useful for understanding the connections between 
solar subsurface properties and magnetic activity in the chromosphere and corona. With the 
improvement of our understanding of acoustic-wave and magnetic-field interactions, and 
also the measurement and inversion techniques, the pipeline codes will be regularly revised. 
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